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Abstract 

Suppose  that  one  knows  a  very  accurate  approximation  a  to  an 
eigenvalue  A  of  a  symmetric  tridiagonal  matrix  T.  A  good  way  to  ap¬ 
proximate  the  eigenvector  x  is  to  discard  an  appropriate  equation,  say 
the  rth,  from  the  system  {T  —  al^x  —  0  and  then  to  solve  the  resulting 
underdetermined  system  in  any  of  several  stable  ways.  However  the 
output  X  can  be  completely  inaccurate  if  r  is  chosen  poorly  and  in  the 
absence  of  a  quick  and  reliable  way  to  choose  r  this  method  has  lain 
neglected  for  over  35  years. 

We  show  how  double  triangular  factorization  (down  and  up),  which 
is  closely  related  to  ‘twisted  factorization’,  gives  us  directly  the  redun¬ 
dancy  of  each  equation  and  so  reveals  the  set  of  good  choices  for  r. 

The  results  extend  to  band  matrices  and  the  applications  go  be¬ 
yond  eigenvector  computation  to  determinant  evaluation  and  solution 
of  well  conditioned  systems. 
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1  Introduction 

The  task  that  started  these  investigations  is  the  computation  of  eigenvectors 
of  a  symmetric  tridiagonal  matrix  (entry  {i,j)  vanishes  if  |i  —  j\  >  1)  once 
the  eigenvalues  are  in  hand.  This  is  not  a  new  problem  and  there  are  good 
programs  available  in  libraries  such  as  LAPACK  and  NAG.  Nevertheless  the 
experts  do  not  consider  the  situation  satisfactory,  see  [9];  the  complexity  of 
the  programs  seems  out  of  proportion  to  the  difficulty  of  the  task  and  the 
adaptation  of  the  current  versions  of  inverse  iteration  to  parallel  mode  is 
frustrating. 

Let  us  briefly  sketch  the  situation.  Given  an  accurate  approximation  a 
to  an  eigenvalue  A  of  an  n  x  n  symmetric  tridiagonal  matrix  T  one  considers 
the  solution  ®  to  the  system  of  equations 

{T  -  al)x  =  h  (1) 

where  h  is  to  be  chosen  wisely.  Since  cr  A  the  best  choice  for  h  is  the 
eigenvector  we  seek  but  this  is  not  an  option.  Next  best  is  to  choose  for  h 
a  column  r  of  the  identity  matrix  I  =  (ei,  62,  • .  • ,  6„).  As  will  become  clear 
in  Section  3  choosing  6  =  is  equivalent  to  omitting  Equation  r  from  the 
system  (1).  The  value  r  =  n  was  proposed  by  Wallace  Givens  in  1954,  see 
[5],  but  no  fixed  value  of  r,  independent  of  a  and  T,  will  do. 

Here  is  a  quotation  from  Wilkinson  concerning  the  computation  of  an 
eigenvector  Uk,  in  Chap.  5,  Section  50,  below  Equation  (50.3)  of  [18]: 

‘Hence  if  the  largest  component  of  Uk  is  the  rth,  then  it  is  the 
rth  equation  which  should  be  omitted  when  computing  Uk.  This 
result  is  instructive  but  not  particularly  useful,  since  we  will  not 
know  a  priori  the  position  of  the  largest  component  of  Ujt.’ 

Ipsen,  in  a  very  readable  survey  attributes  the  idea  of  omitting  one  equa¬ 
tion  of  the  system  to  Wilkinson,  see  Section  7  of  [9],  but  we  suspect  that 
this  method  was  routinely  taught  in  mathematics  classes  before  Wilkinson 
was  born,  see  [8],  [2],  and  [14].  He  was  born  in  1919  and  [8]  was  published 
in  1921. 

Wilkinson  abandoned  the  hunt  for  a  good  value  of  r  and  used  b  =  PLe 
where  T  -  al  =  PLU  denotes  triangular  factorization  with  partial  pivoting 
and  e  =  ±e,-,  see  [17].  However  even  this  choice  fails  if  some  eigenvalues 
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axe  equal  to  working  accuracy  and  he  resorted  to  ‘tweaking’  the  computed 
eigenvalues  in  such  cases. 

In  private  communication  to  one  of  us  Wilkinson  declared  that  he  would 
prefer  6  =  to  6  =  PLe  if  only  he  knew  a  quick,  reliable  way  to  choose  r 
so  that  the  rth  entry  of  the  wanted  vector  is  above  average,  not  necessarily 
the  greatest. 

The  current  LAPACK  codes,  see  [4],  do  not  use  Wilkinson’s  choice;  in¬ 
stead  b  is  chosen  ‘at  random’  from  an  appropriate  distribution  but  this  makes 
it  difficult  to  obtain  orthogonal  eigenvectors  for  close  eigenvalues.  The  case 
for  this  approach  is  made  in  [10]. 

In  this  paper  we  present  a  new  way  to  choose  r  that  depends  strongly  on  T 
and  a.  However  it  is  not  free;  the  cost  is  essentially  n  extra  divisions.  It  turns 
out  that  our  method  produces  further  information,  beyond  the  right  value 
of  r,  that  helps  us  avoid  the  computation  of  completely  negligible  entries  in 
the  wanted  eigenvector.  In  this  way  the  overhead  for  finding  the  right  value 
of  r  pays  for  itself  as  n  becomes  large.  See  Figure  4  after  reading  Section  3. 

Our  method  uses  two  complete  triangular  factorizations,  one  starts  from 
the  top  and  the  other  from  the  bottom.  This  idea,  of  itself,  is  not  new  and 
forms  the  basis  of  ‘twisted  LDU\  What  has  not  been  noticed  before  is  that 
by  combining  both  sets  of  ‘pivots’  one  finds  the  redundancy  measure  of  each 
row.  Then  one  is  in  a  good  position  to  choose  r.  Twisted  factorization, 
in  contrast,  stops  the  eliminations  when  they  meet  at  some  predetermined 
interior  row.  By  completing  the  up  and  down  factorizations  at  a  total  cost 
of  2n  divisions  we  have  full  information  on  all  possible  twisted  factorizations 
each  of  which  costs  n  divisions.  A  few  historical  remarks  on  twisted  LDU 
are  given  at  the  end  of  Section  3. 

Section  2  discusses  the  ‘obvious’  solution  to  the  problem  and  shows  its 
shortcomings.  The  new  method  is  implicit  in  Theorem  1  which  is  established 
in  Section  3  along  with  Theorem  2  which  presents  accurate  ways  to  compute 
the  determinant.  Section  4  shows  how  the  quantities  introduced  in  Theo¬ 
rem  1  reveal  the  ‘envelope’  of  an  eigenvector  when  the  tridiagonal  is  normal. 
Section  5  extends  the  results  to  cover  breakdown  in  triangular  factorization 
and  zero  entries  in  eigenvectors.  Section  6  extends  the  results  of  Section  3  to 
block  tridiagonal  matrices.  Application  of  these  ideas  and  error  analysis  will 
be  given  elsewhere. 

The  reader  is  expected  to  know  the  LDU  theorem  concerning  existence 
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and  uniqueness  of  triangular  factorization  and  the  expressions  for  the  pivots, 
as  the  diagonal  entries  of  D  are  often  called.  In  this  representation  both  L 
and  U  have  I’s  on  the  diagonal.  In  practice,  when  division  is  slow,  people 
often  use  [LD)D  ^  {DU)  instead  of  LDU  but  the  distinction  is  not  important 
in  this  paper. 

The  main  notational  issue  is  the  representation  of  submatrices. 

In  MATLAB  notation  the  submatrix  of  M  in  rows  i  through  j  and  columns  k 
through  1  is  given  by  M{i  :  j,k  :  1).  This  is  clear  but  sometimes  too  obtrusive. 
We  use  to  denote  the  principal  submatrix  M{i  :  j,i  :  j).  For  column 
vectors  we  prefer  simple  lower  case  Latin  letters  aj,  y, . . .  in  bold  face  type, 
with  entries  x(l),  x{2), . . . ,  a:(n).  For  subvectors  we  use  either  x{i  :  j)  or  x''U 
Finally  we  try  to  use  lower  case  Greek  letters  a,  /3, . . .  for  scalars  although 
matrix  entries  will  be  written  as  M{iJ)  or  Mij. 

One  notational  innovation  is  to  use  +  to  indicate  a  process  taking  rows 
in  increasing  order  and  —  to  indicate  the  process  going  in  decreasing  order, 
e.  g.  LDU  is  written  as  L+D+U^  while  UDL  is  written  as  U-D-L^. 

As  usual  ||i;||  =  ||r||2  =  v^,  while  \\v\\^  =  max,-  |u(i)|.  The  dimension, 
or  order,  of  Ci  or  of  any  column  of  the  identity  matrix  I  =  (cj,  63, ... ,  e„)  is 
given  by  the  context. 

Theorem  1  was  presented  at  the  SIAM  conference  on  Parallel  Processing 
in  San  Francisco,  California  in  February  1995  by  two  of  us  (KVF  and  BNP). 
Previous  work  that  used  a  different  method  to  compute  the  of  Theorem  1, 
a  technique  more  prone  to  overflow,  was  presented  (  by  KVF  and  BNP)  at 
the  Householder  XII  conference  at  Lake  Arrowhead,  California  in  June  1993, 
and  the  SIAM  Applied  Linear  Algebra  meeting  in  Snowbird,  Utah  in  June 
1994. 

2  A  Classical  Analysis 

In  case  a  pure  mathematician  should,  by  chance,  read  this  material  it  seems 
wise  to  begin  by  explaining  that  the  problem  discussed  here  is  not  as  trivial 
as  it  may  appear  at  first.  It  is  the  computer’s  limited  precision  that  causes 
the  diflSculties. 

Anyone  who  has  mastered  an  introductory  course  in  matrix  theory  and 
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who  h^ls  absorbed  the  significance  of  the  tridiagonal  form  J  (with  nonzero 
values  adjacent  to  the  diagonal)  might  reason  as  follows. 

Lemma  1  An  eigenvector  of  an  unreduced  tridiagonal  matrix  J  cannot  have 
a  0  in  the  first  or  last  component. 

Proof.  Consider  the  equation  for  an  eigenvector  x  {jfi  0)  associated  with  an 
eigenvalue  A, 

( J  -  \l)x  =  0.  (2) 

Suppose  that  a:(l)  =  0.  Then  the  first  equation  in  (2)  dictates  that  x{2)  = 
(A  — Jii)a;(l)/Ji2  =  0  as  well,  since  Jn  ^  0.  Now  the  second  equation  dictates 
that  a:(3)  is  a  linear  combination  of  x(l)  and  x{2)  and  also  vanishes.  Pro¬ 
ceeding  with  the  remaining  equations,  in  order,  it  appears  that  every  entry 
of  X  must  vanish  in  contradiction  to  the  assumption  that  x  is  an  eigenvector. 
So  the  assumption  that  a:(l)  =  0  is  not  tenable.  By  similar  reasoning  but 
taking  the  equations  in  reverse  order  it  is  untenable  that  x(n)  =0.  □ 

The  preceding  argument  also  shows  one  way  to  compute  an  eigenvector  of 
J.  It  is  valid  to  set  x(l)  =  1  and  to  use  the  first  equation  of  (2)  to  determine 
x(2),  and  the  second  to  determine  x(3),  using  x(l)  and  x{2).  Proceeding 
as  before  the  rth  equation  may  be  used  to  determine  x(r  +  1)  and  thus  x 
may  be  obtained  without  actually  making  use  of  the  nth  equation  which, 
says  the  mathematician,  will  be  satisfied  automatically  since  the  system  (2) 
is  singular. 

It  would  be  equally  valid  to  begin  with  x(n)  =  1  and  to  take  the  equations 
in  reverse  order  to  compute  x{n  —  1), ... ,  x(2),  x(l)  in  turn  without  using  the 
first  equation  in  (2).  When  normalized  in  the  same  way  x  and  x  will  yield 
the  same  eigenvector.  Note  that  the  problem  has  been  solved  without  the 
bother  of  computing  a  triangular  factorization. 

The  proof  of  Lemma  1  actually  shows  a  little  more  than  was  claimed. 
For  an  upper  Hessenberg  matrix  {{i,j)  entry  vanishes  if  t  >  j  -|-  1)  that  is 
unreduced  (entries  {i  -f  l,i)  do  not  vanish)  x{n)  cannot  vanish  and  for  an 
unreduced  lower  Hessenberg  matrix  x(l)  cannot  vanish. 

The  method  described  above  was  proposed  by  W.  Givens  in  1954,  see  [5]. 
It  often  gives  good  results  when  realized  on  a  computer  but,  at  other  times, 
delivers  vectors  pointing  in  completely  wrong  directions. 
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The  preceding  analysis  is  valid  in  exact  arithmetic  but  is  inapplicable  to 
computer  work  for  the  following  reasons.  First,  it  is  rare  that  an  eigenvalue 
of  a  tridiagonal  (or  any  other)  matrix  is  representable  in  limited  precision. 
Consequently  the  systems  such  as  (2)  that  are  to  be  solved  are  not  singular 
and,  in  (2),  the  unused  equation  will  not  be  satisfied  automatically  even  if  the 
solutions  of  the  other  equations,  in  turn,  were  obtained  exactly.  The  second 
weakness  is  that,  in  a  computer,  the  sequence  1,  x(2), . . . ,  can  overflow.  This 
is  a  possibility  that  pure  mathematicians  do  not  have  to  worry  about. 

It  turns  out  that,  for  isolated  eigenvalues,  Givens’  method  gives  an  excel¬ 
lent  approximate  eigenvector  whenever  the  first  or  last  entry  of  the  wanted 
eigenvector  is  above  average  in  magnitude.  Conversely  it  gives  disastrous  re¬ 
sults  when  those  extreme  entries  are  tiny.  Wilkinson  gives  a  striking  example 
in  Section  52,  Chap.  5  of  [18]. 

The  purpose  of  this  section  was  to  show  that  the  ‘obvious’  method  for 
computing  eigenvectors  is  not  adequate  for  finite  precision  arithmetic. 

3  Diagonal  of  the  Inverse 

In  basic  courses  in  matrix  theory  one  is  taught  to  solve  a  system  of  equations 
by  computing  a  row  echelon  form.  If  the  system  is  singular  at  least  one 
row  of  the  echelon  form  vanishes  and  the  corresponding  row  of  the  original 
system  is  redundant.  The  homogeneous  system  is  solved  by  assigning  any 
values  to  the  ‘free’  variables  and  backsolving  for  the  rest  of  them.  In  general 
a  discarded  row  is  not  unique;  it  need  only  be  a  linear  combination  of  the 
remaining  ones. 

In  practice  our  system  is  nearly,  but  not  quite,  singular  and  a  natural 
modification  of  the  standard  procedure  is  to  seek  a  row  that  is  most  nearly 
redundant  and  then  ignore  it  while  determining  a  solution  x  to  the  remaining 
homogeneous  system.  This  solution  x  will  not  satisfy  the  omitted  (rth) 
equation.  In  other  words,  faced  with  the  fact  that  Mx  =  0  admits  only  the 
trivial  solution  one  finds  a  suitable  r  and  solves,  instead,  Mx  =  e^Sr  where 
is  the  ‘defect’  or  residual  of  the  rth  equation. 

This  is  what  is  meant  by  ‘omitting  the  rth  equation’. 

In  general  it  is  difficult  to  find  r  and  to  solve  the  reduced  homogeneous 
system.  Fortunately  when  M  is  tridiagonal  the  omission  of  row  r  splits  the 
system  into  two  separate  parts.  For  a  modest  cost  the  residual  for  every 
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choice  of  can  be  computed  and  that  gives  an  excellent  basis  for  choosing 
the  right  r. 


Theorem  1  (Double  Factorization)  Let  J  he  a  tridiagonal  nxn  complex 
matrix  that  permits  triangular  factorization  in  both  increasing  and  decreasing 
order  of  rows: 


L+D+U+  =  J=U.D.L.. 

(3) 

For  each  k^ 

1  <  k  <  n,  define  ■jk  o.'nd  by 

=  6*7*,  z^'‘\k)  =  1. 

(4) 

Then 

7*  =  D^{k)  +  D-{k)  —  Jkk- 

(5) 

Proof  In  what  follows  MATLAB  notation  will  be  used  for  submatrices  that 
are  not  square  and  a  more  condensed  representation  otherwise.  In  addition, 
if  terms  that  involve  out  of  range  indices  are  dropped  then  the  analysis  that 
follows  covers  the  extreme  cases  k  =  1  and  k  =  n  a^s  well.  For  brevity  write 
for 

Omit  the  kth  equation  from  (4)  and  what  remains  is  two  homogeneous 
systems.  Next  use  the  appropriate  triangular  factorization  (3)  to  write  these 


systems  as 

1  .  jt)z(l  I  k)  =  0,  (6) 

^fc+i:n^*+i:n^_ (jfc  +  1  ;  :  n)z(ifc  I  n)  =  0.  (7) 

By  the  assumption  that  the  LDU  and  U DL  factorizations  exist  the  matrices 

must  be  invertible.  Premultiply  (6)  and  (7)  by 
the  appropriate  inverses  to  find 

7/+(l  :  A:-l,l  :  A:)z(l  :  fc)  =  0,  (8) 

L-{k  +  1  :  n,k  :  n)z{k  :  n)  =  0.  (9) 

The  last  equation  in  (8)  shows  that 

1  •  - 1)  +  Uj^-i,k4k)  =  0.  (10) 

The  first  equation  in  (9)  shows  that 

1  +  1)  =  (11) 
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Recall  that  z[k^  —  1  and  substitute,  from  (10)  and  (11),  the  values  for  z(ifc — 1) 
and  z{k  +  1)  into  the  fcth  equation  of  (4)  to  find,  for  Jb  =  2, . . . ,  n  —  1, 

"b  '^k,k  ~  'Jk,k-\-lLf._^^  f^  (12) 

=  {Jkk  —  Jk,k-\Uk_-^,^)  —  Jk,k  +  {Jkk  —  Jk,k+\Llj^-^k) 

=  Dj^{k)- Jkk-\- D.{k),  (13) 

as  claimed.  For  A:  =  1  note  that  I>+(1)  =  Ji,i  and  71  =  D.{1).  For  k  =  n 
note  that  D-(n)  =  J„n  and  7„  =  D+(n).  Thus  (5)  holds  for  A:  =  1  an  A:  =  n 
as  well  as  for  A:  =  2, . . . ,  n  —  1.  □ 

Corollary  1  Let  J  satisfy  the  Hypotheses  of  Theoretn  1.  Either  J  is  singular 
and  then  D+{n)  =  T)_(l)  =  7^  =  7j  =  0  and  both  and  are  in  J’s 
null  space,  or 

diag{J~'^)~'^  +  diag{J)  =  D+  +  D-  .  (14) 

Proof  By  the  assumption  of  (3)  the  only  D  values  that  can  vanish  are  D+{n) 
and  D_(l).  Also  T>+(1)  =  Ju  and  D-{n)  =  so  that,  when  J  is  singular 

7n  =  0  -  +  D.{n)  =  0,  7i  =  £>+(1)  -  +  Q  =  0, 

and  so  Jz^^l  =  Jz^"!  =  0. 

If  J  is  invertible  then  7*  must  be  nonzero  for  all  A:  =  1, ...  ,n  (to  avoid 
giving  a  nontrivial  solution  to  Jx  =  0)  and  multiplication  of  (4)  by  7^ 
yields 

Ik^  =  Ik^z^^Kk)  =  =  elJ-^ek.  (15) 

Thus 

'yk  ^  7— IN  ’  ^  1,  .  .  .  ,  71. 

W  jkk 

□ 

Equation  (14)  is  a  striking  property  of  invertible  tridiagonals  and  gave  us 
the  title  for  this  section. 

In  applications  it  is  useful  to  have  several  different  expressions  for  in 
addition  to  (13). 
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Corollary  2  With  the  notation  of  the  Theorem  1,  for  1  <  A:  <  n, 

D+{k)-Jk,k  +  D.{k), 

+  Jk,k  -  Jk-\-\,kU]f^k+\-> 

~Jk,k-lU^_l  li.  +  Jk,k  ~  ‘^fc,fc+l-^fc+l,fc» 

-^+(^)  ~  ^k,k+l'^k+l,k, 

~^k,k-l'^k-l,k  +  D-{k). 

For  k  =  1  and  k  =  n  omit  terms  with  invalid  indices. 

Proof.  The  first  and  third  expression  are  just  (5)  and  (12).  The  others  come 
from  rewriting  (12)  as 

Ik  =  —Jk,k-iJk-i,k/F+{k  —  1)  +  Jk,k  —  Jk,k+iJk+i,k/D+{k  +  1)  (16) 

and  using  Jk,k+i  =  UlMD.(k  +  1)  =  D^(k)Uli,^,  etc.  and  D.(k)  = 
Jkk  —  Jk,k+iJk+i,k/D-{k  +  1),  etc.  □ 

When  J  is  nearly  singular  one  is  most  interested  in  the  values  of  k  that 
yield  minimal  |7fc[  values. 

The  middle  formula  in  Corollary  2  is  of  most  interest  to  us  because  of  the 
following  result  which  shows  that  no  divisions  are  needed  to  find  once  k 
is  known.  . 

Corollary  3  With  the  notation  of  the  Theorem  1  and  a  given  value  of  k,  so 
that  z  =  Jz  =  Ck'ik,  then 

z{j)  =  -C^i!'j+i^(i  +  1),  i  =  A;  - 

•2^(*)  f ))  i  —  k  1,  ,  .  .  .  ,  Tl. 

Proof  These  equations  are  (8)  and  (9)  in  expanded  form.  □ 

Another  reward  for  computing  both  factorizations  is  a  wide  choice  of 
expressions  for  det  J. 


10 


Theorem  2  Assume  the  hypothesis  of  Theorem  1.  Then  for  =  1, . . .  ,n, 


r  £>+(1)  •  •  •  D^{k  -  l)^kD.{k  +  1)  •  ■  •  D_(n) 


det  J  =  < 


D^{k-2)det 


D+{k  —  1)  Jk-i,k 
Jk,k-i  D.ik) 


D-{k  +  !)•••  D-{n) 


Tk  _  D+{k) 

Tk+i  D-{k  +  1) 


Proof.  Apply  Cramer’s  rule  for  the  ^th  entry  of  where  =  Ckjk- 

The  numerator  is  a  determinant  whose  kth  column  is  ek'fk-  Expand  it  by 
column  k  to  find 


1  =  z(>^){k)  =  jk  det  det  j'‘+^-'^/det  J. 

Since  J  =  L^D^Ujf.  =  it  follows  that 

det  =  T>+(1)  • .  •  D+{k  -  1),  det  =  D.{k  +  l)---  D.{n). 

The  second  expression  comes  from  the  twisted  factorization  of  J: 


■  0 

n 

■  0 

0  C/^=” 

LJ 

£>i+l:n 

0  X*=" 

where 

r£)+(A:-l)  Jk-i,k 

[  Jk,k-i  D.{k)  • 

From  the  first  expression  for  det  J  it  follows  that 

')kD-{k  +  1)  =  7;b+i£>+(A:), 

which  gives  the  ratio  of  consecutive  7’s.  □ 

When  there  is  severe  cancellation  in  computing  7^  from  any  of  the  formu¬ 
lae  in  Corollary  2  then  it  may  be  possible  to  take  extra  care  in  the  evaluation 
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of  det  □  and  win  back  a  few  bits  of  precision.  If  warranted  the  idea  may  be 
taken  further  to  use 


det  tridiag 


D+{k  —  1)  Jk-i,k  0 
Jk,k-1  Jk,k  Jk,k+l 

0  Jk+l,k  D-(k  +  1) 


for  the  sensitive  part  of  the  computation.  These  details  are  of  great  practical 
importance  when  J  is  close  to  singular  as  occurs  in  iterative  methods  for 
finding  eigenvalues. 


Remark  1  An  attentive  reader  may  be  puzzled  that  Corollary  3  cannot 
generate  an  isolated  0  entry  in  If  z{j)  vanishes  (because  J{j,j  + 1)  =  0) 
then  all  entries  z{l),  I  <  j  <  k,  must  vanish  too.  This  is  appropriate  since  the 
matrix  is  reduced  when  Jjj+i  =  0.  Yet  there  exist  eigenvectors  of  unreduced 
tridiagonals  with  isolated  entries  that  vanish;  such  entries  are  the  nodes  of 
the  eigenvector.  The  explanation  is  that  the  hypothesis  (13)  does  indeed  rule 
out  isolated  zero  entries.  Section  5  extends  the  results  of  Theorem  1  to  cover 
these  important  cases.  There  we  see  that  hypothesis  (3)  in  theorem  1  is  not 
essential. 


Remark  2  Corollaries  2  and  3  show  that  we  need  only  retain  the  nontrivial 
entries  of  L-  and  U+  in  the  factorization  process  in  order  to  obtain  all  the 
7- values,  and  for  any  given  r,  to  solve  for  with  no  more  divisions. 

Remark  3  For  large  n  there  will  be  many  products  in  the  calculation  of 
z{\)  or  z{n)  for  a  given  r.  In  general  one  is  concerned  about  possible  overflow 
but  here  r  is  selected  so  that  z{r)  should  be  a  maximal,  or  nearly  maximal, 
entry  of  z  and  so  no  overflow  can  occur.  Underflow,  if  it  occurs,  is  harmless 
here  and  should  be  flushed  to  0. 

Remark  4  From  (4)  we  see  that  the  FP  vector  is  annihilated  exactly 
by  J  —  6^7,6*,  a  rank  one  perturbation. 

Let  us  summarize  this  section  in  the  form  of  a  high  level  algorithm. 


Null  Vector  Approximation. 
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Input:  J,  nearly  singular,  n  x  n,  tridiagonal 

1.  Factor  J  down  and  up  to  compute  -y  =  (71, . . .  ,7„)  using  a  convenient 
formula  from  Corollary  2. 

2.  Find  min,- 17,-|,  perturbing  zero  values  appropriately.  Compute  the  index 

set  C  :=  :  |7_,-|  <  fi  min,-  [7^!  }.  Choose  fx  so  that  1  <  fx  <2. 

3.  Select  a  suitable  r  from  C  and  compute  using  Corollary  3.  It  is 
also  possible  to  take  an  appropriate  linear  combination  of  two  or  more 

j  €  C. 

Output:  r  and  z^’’^ 

We  do  not  claim  that  z^’'),  which  we  call  the  FP  vector,  is  always  an 
adequate  approximation  to  a  null  vector  of  J .  However  it  is  always  useful.  We 
do  not  discuss  the  calculation  of  orthogonal  vectors  for  clustered  eigenvalues 
here,  see  [13].  It  seems  to  be  wise,  when  eigenvalues  are  real,  to  forbid  the 
same  r  to  adjacent  eigenvalues. 

There  are  several  possibilities  for  replacing  zero  values  of  'fk  in  Step  2. 
Theorem  2  gives  one  way  and  Remark  4  suggests  that  we  set  7*  =  macheps  ■ 

Jkk- 

Estimates  of  the  accuracy  of  as  an  approximate  null  vector  may  be 
given  in  different  contexts:  the  nonsymmetric  case,  the  real  symmetric  case, 
the  general  linear  eigenvalue  problem,  and  the  bidiagonal  singular  vector 
case.  We  do  not  wish  to  submerge  the  ideas  in  this  paper  with  such  results 
but  see  the  corollary  to  Theorem  3  in  Section  4. 

As  mentioned  in  the  introduction  twisted  LDU  starts  elimination  from 
the  top  and  the  bottom  and  stops  at  some  selected  interior  row  k.  Henrici 
(1963)  used  twisted  LDU  implicitly  in  deriving  optimal  Gersgorin  bounds,  in 
[7],  for  unreduced  real  tridiagonals  with  some  complex  eigenvalues.  D.  Ker¬ 
shaw  (1970),  in  [11],  obtained  nice  bounds  on  ( J~^)rr/ Jrr  using  twisted  LDU. 
Babuska  (1972)  wanted  a  specific  entry  in  J~^h  for  any  h  and  his  formula 
(5.29)  on  page  62  of  [1]  is  one  instance  of  our  formula  (5).  Fischer  et.  al. 
(1974),  in  [6],  discuss  a  twisted  Toeplitz  factorization  of  Buneman  and  at¬ 
tribute  the  adjective  ‘twisted’  to  Strang  [15].  Dongarra  et.  al.  (1979),  in 
the  LINPACK  codes  use  twisted  LDU  meeting  in  the  middle  for  improved 
efficiency  and  the  practice  has  been  taken  up  in  parallel  computation,  see 
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[16].  They  call  it  the  BABE  algorithm  (Begin,  or  Burn,  At  Both  Ends),  see 
[3]. 


4  The  Eigenvector  Connection 

Suppose  that  J  is  a  normal  matrix  as  well  as  tridiagonal; 

j  =  VAV\  =  V*, 


and 

A  =diag{Xi,...,Xn),  V  =  (vi, . . . , 

Theorem  3  Let  J  =  J  —  al  satisfy  the  hypotheses  of  Theorem  1.  Then 
=  7fc(cr)^  for  each  k,  and  as  a  — ►  where  Xj  is  an  isolated  eigenvalue 

of  J y  then 

tm 

for  all  k,  1  <  k  <  n,  such  that  Vj{k)  ^  0. 

Proof.  From  (15),  for  each  k 


=  efc(J-<T/)  ^Bk 

=  elV{K-cI)-^V^ek 

h 

=  ih  -  ^  iv.(^)p  I  • 

Sum  for  A:  —  1, . . . ,  n  and  use  the  orthogonality  of  V  to  find 


n 


(A,  -  <7)-^ 

(Xj  -  a)"^ 
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Hence,  as  a  — y  Aj,  an  isolated  eigenvalue, 


E"  /v-l 
m=l  Im 

provided  that  >  0. 


—  a 

(j 


i  +  E 


Aj  - 
\i  -  c , 


-1 


□ 


Since 

E|..,(i)p  =  i-M*)p 

‘Yi 

convergence  is  more  rapid  the  larger  is  |ui(fc)|.  For  a  given  a,  much  closer 
to  Aj  than  to  any  other  A,-,  the  vector  with  components 
A:  =  1, . . . ,  n,  provides  a  useful  envelope  of  the  eigenvector  Vj.  It  is  pleasing 
that  there  is  no  requirement  that  the  7fc  be  real;  the  appropriate  quotients 
will  be  positive. 

The  FP  vector  z  defined  by  (4)  is  an  alternative  approximation  to  the 
eigenvector  Vj  of  J.  Its  quality  is  indicated  in  the  next  result. 

Corollary  4  With  the  notation  of  Theorem  3  let  Jz  =  (J  —  al)z  =  6^7*. 
Then 

(a)  Rayleigh  quotient{z)  =  a  +  7ife/i|*|P, 

(b)  |Aj  -a-  7fc/||z||^|  <  |7;t|/||2||  •  min{l,  \lk\l{\\z\\gap)] , 

where  gap  =  min{|Aj+i  -  a\,  \a  -  Aj_i  |}, 

(c)  sin Z(vj,z)  <  \lk\l{\\z\\gap). 

Proof  z*Jzlz*z  =  (T-f-z*efc7fc/||z||2,  gives  (a).  Also  \\Jz\\l\\z\\  =  |7/t|/||«|| 
and  standard  gap  theorems,  see  Chapter  11  in  [12],  give  (b)  and  (c).  □ 

Unfortunately,  in  finite  precision  arithmetic,  we  cannot  let  a  — y  Aj  and 
so,  in  practice,  a  7j  may  vanish  in  the  same  way  that  a  divided  difference 
may  vanish  even  when  the  desired  derivative  does  not.  However  we  do  not 
need  more  than  a  few  bits  of  accuracy  in  7j,  it  is  its  exponent  that  counts. 

Figure  1  shows  the  7-vector  for  different  values  of  a  and  the  true  profile 
of  a  simple  eigenvector,  all  on  a  log  scale. 


Theorem  4  An  unreduced  tridiagonal  matrix  is  normal  if,  and  only  if,  it  is 
a  translate  of  a  Hermitian  or  skew~Hermitian  matrix,  possibly  multiplied  by 
a  rotation  factor  e'^,  i^  —  —1. 

The  proof  is  left  to  the  interested  reader. 


5  Zero  Pivots 

Triangular  factorization  is  said  to  fail,  or  not  exist,  if  a  zero  ‘pivot’,  D^{j) 
or  D.{j)  is  encountered  prematurely.  The  last  pivot  is  allowed  to  vanish 
because  it  does  not  occur  as  a  denominator  in  the  computation. 

One  of  the  attractions  of  an  unreduced  tridiagonal  matrix  is  that  the 
damage  done  by  a  zero  pivot  is  localized.  Indeed,  if  oo  is  added  to  the  number 
system  then  triangular  factorization  cannot  break  down  and  the  algorithm 
always  maps  J  into  unique  triplets  L,  D,  U.  There  is  no  need  to  spoil  the 
inner  loop  with  tests.  It  is  no  longer  true  that  LDU  =  J  but  equality  does 
hold  for  all  entries  except  for  those  at  or  adjacent  to  any  infinite  pivot. 

It  is  possible  to  work  with  signed  oo  (affine  geometry)  or  unsigned  oo  (the 
complex  plane)  and  it  will  be  easiest  for  our  purposes  to  use  the  unsigned 
oo.  Thus  +1/0  =  —1/0  =  oo. 

If  we  allowed  off  diagonal  entries  to  vanish,  in  which  case  J  is  said  to  be 
reduced,  then  we  might  encounter 

T(^  +  1,  A:)  =  J(A:  +  1,  fc)/Z)(ifc)  =  0/0 

and  that  would  be  a  genuine  breakdown. 

Let  us  examine  what  happens  when  D{k  —  1)  =  0.  In  turn 

L{k,k-\)  =  J{k,k-l)lD{k-l)  =  oo, 

•  U{k-\,k)  =  J{k-\,k)ID{k-l)  =  oo, 

JD{k)  =  J(^k,k')  L{k,k  —  1)T(A:  —  l,A:)  =  oo, 

L(^  +  1,A;)  =  j(^kAl,k)ID{k)  =  Q, 

•  U{k,k^\)  =  J{k,k  +  l)ID{k)  =  Q, 

D{k  +  1)  =  ^(A:  +  l,A;  +  l)-Z(A:+l,)t)J(jfc,ib  +  l) 

=  J(A;  +  1,  A:  +  1). 
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Unless  J(Jk  +  l,/:  +  l)  =  0  the  factorization  proceeds  normally  until  the  next 
zero  pivot  is  encountered.  We  have  placed  an  •  against  entries  that  are  not 
computed  when  a  simple  LU  factorization  is  used.  Here  U  =  DU  in  the 
finite  case. 

When  the  product  LDU  is  formed  in  the  case  given  above  then  various 
strange  expressions  such  0  •  oo  and  cx)  +  oo  arise  and  we  designate  them 
by  NaN  (Not  a  Number).  We  discover  that  LDU  =  J  except  in  row  and 
column  k.  Note  that  D{k)  =  oc. 

It  is  important  to  our  later  results  to  show  that  when  J  is  singular  then 

D^{k)  =  oo  if,  and  only  if,  D^{k)  =  oo, 

where  the  notation  follows  Section  3. 

It  turns  out  that  infinite  pivots  correspond  to  zero  entries  in  eigenvectors 
and  so  have  a  legitimate  role  in  the  theory. 

Theorem  5  Let  J  be  n  x  rij  tridiagonaly  unreducedy  and  singular.  For  each 
ky  \  <  k  <  Uy  is  singular  ify  and  only  ify  is  singular.  They  are 

singular  ify  and  only  ify  z{k)  =  0  whenever  Jz  =  0. 

Proof,  Write 


and  partition  Jz  =  0  conformably.  Thus 

+  Jk-i,kz{k)ek.i  =  0,  (17) 

eiJk+i,kz{k)  +  =  0,  (18) 

and  «+(!)  7^  0,  z+(n)  ^  0  by  Lemma  1  in  Section  2. 

If  z{k)  =  0  then  (17)  shows  that  z+(^  0)  is  in  null  space  and  (18) 

shows  that  z-{^  0)  is  in  null  space.  So  both  matrices  are  singular. 

Now  consider  the  converse,  z{k)  7^  0.  Since  J  is  unreduced  rank{J)  = 
n  —  1  and  its  null  space  is  one  dimensional.  So  the  system 

Jz  =  0,  z{k)  =  1, 

has  a  unique  solution.  Thus  both  (17)  and  (18)  are  inhomogeneous  equations 
with  unique  solutions.  Thus  and  are  invertible.  □ 


17 


Corollary  5  Let  J  be  n  y.  n,  tridiagonal,  unreduced,  and  singular.  Let  the 
triangular  factorization  algorithm  applied  to  J  in  both  increasing  and  decreas¬ 
ing  order  of  rows  yield  unique  matrices  C/+,  [/_,  D_,  and  X_.  Then, 

for  j  =  1,2,..., n, 

I>+(j)  =  oo  iff  Z)_(»  =  oo. 

Proof. 

£>+(;)  =  oo  -  1)  =  0 

4=^  singular 

ji+i:n  singular  (by  Theorem  5) 

D-{j  +  1)  =  0 

D-{j)  =  oo. 

□ 

In  Theorem  1  of  Section  3  the  value  of  7*  was  fixed  by  the  condition 
z{}t)  =  1  imposed  on  the  solution  oi  Jz  =  6^7*-  When  J  is  singular  there  is 
a  nonzero  solution  to  Jz  =  0  and  the  attempted  normalization  z{k)  =  1  is 
valid,  even  if  not  wise,  in  all  cases  except  when  z{k)  =  0. 

An  appropriate  signal  that  an  infeasible  normalization  has  been  imposed 
is  that  7fc  =  NaN  (Not  a  Number)  and  that  is  precisely  what  the  formulae 
in  Corollary  2  deliver  whenever  J  is  singular  and  D+{k)  =  D-{k)  =  00.  In 
these  cases,  in  exact  arithmetic,  D+{n)  =  0  and  7„  =  0  as  well  as  £>-(!)  =  0 
and  7i  =  0.  Thus  in  the  search  for  a  minimum  value  of  |7_,|  indices  j  that 
have  7j  =  NaN  will  never  be  selected. 

The  good  news  is  that  by  computing  all  the  {7^}  it  is  known  in  advance 
whether  or  not  the  z^*')  has  a  zero  entry.  In  the  generic  case,  with  no  zeros, 
the  algorithm  given  in  Corollary  1  in  Section  3  may  be  used  free  of  any  tests 
for  invalid  operations.  In  the  exceptional  case  the  following  procedure  may 
be  used. 
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I  -  1),  ■2(*  -  1)  7^  0,  ) 

\  -(Ji-t,i-2z(i  -  2)fJi-i,i),  otherwise  J 


i  =  fc  +  1, . . . ,  n 


This  algorithm  will  not  touch  any  infinite  values  in  [/+  or  L_. 

When  J  is  not  singular  Theorem  1  continues  to  hold,  if  oo  is  allowed,  in 
the  following  sense. 


Corollary  6  (of  Theorem  2)  If  J  is  unreduced,  tridiagonal,  and  oo  is  rep¬ 
resented  then 


D+{k  -  1)  =  0  implies  (J  ^)kk  -0  (7*  =  00), 

D.{k)  =  0  implies  =  0  (7fc_i  =  00). 

Proof,  Use  the  twisted  factorization  in  the  proof  of  Theorem  2  that  introduces 
the  2x2  matrix 

p-j  _  D^(k  “  1)  Jk-i,k 

[  Jk^k^i  D^{k)  _  * 

Invert  J  and  observe  that  there  is  a  simple  expression  for  the  {k^k)  and 
(fc  —  1,  —  1)  entries: 

=  (□-')2.2  =  D+{k  -  l)/det  □  . 

If  J  is  unreduced  and  D^{k  1)  =  0  then  dei  □  =  —Jk,k-iJk-\,k  ^  0.  This 
establishes  the  first  aissertion.  Similarly 

{J-%-r,k-i  =  D.{k)ldetU. 

□ 


Figures  2  and  3  show  a  striking  instance  of  Theorem  5  for  the  matrix 
W21  and  the  pair  of  eigenvalues  close  to  6.  Each  horizontal  line  of  the  figure 
corresponds  to  one  value  of  k]  eigenvalues  of  are  marked  by  +  and 

eigenvalues  of  are  marked  by  o.  Theorem  5  implies  that  if  an  eigen¬ 

vector  has  a  zero  entry  in  position  k  then  a  o  and  a  -1-  must  coincide  on  the 
eigenvalue  in  line  k.  Indeed,  in  Figure  3  (an  enlarged  picture  of  Figure  2  near 
6),  when  fc  =  11  this  is  precisely  what  happens.  For  neighboring  values  of  k 
the  Ritz  values  are  not  particularly  close  to  eigenvalues  and  after  it  =  1 1  the 
o  is  replaced  by  a  -f-  in  the  interval  (A12,  A13).  If  v  is  a  normalized  eigenvector 
with  eigenvalue  A  then  v{ky  is  proportional  to  the  product  of  the  distances 
of  A  from  the  -1-  and  o  points  on  line  k. 
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6  Block  Tridiagonals 

If  an  arithmetic  system  lacks  the  symbol  oo  it  is  possible  to  extend  Theorem  1 
by  using  blocks  in  the  LDU  factorization.  If  J  is  unreduced  then  there  always 
exists  a  factorization 


L+D+U+  =J=U-D.L- 

if  the  D's  are  allowed  to  have  2x2  and  1x1  blocks  along  diagonal,  no  larger 
blocks  are  needed.  However  Theorem  1  extends  beyond  this  case  to  band 
matrices  and,  to  any  block  tridiagonal  matrix.  Thus  Z)+  and  Z)_  are  direct 
sums  of  square  blocks;  T+  and  17+  are  conformable  with  £>+,  L_  and  U-  are 
conformable  with  jD_. 

Theorem  6  Let  J  permit  block  triangular  factorization  in  both  increasing 
and  decreasing  order  of  indices 


L+D+U+  =  J=  U.D.L. . 


There  is  no  requirement  that  the  block  structures  of  D+  and  D-  be  con¬ 
formable.  However  for  any  corresponding  blocks  k  and  I  such  that  D+{k) 
and  D-{1)  are  conformable  and  m  by  m,  define  the  m  x  m  matrix  F  by 
equations 


Z+  N 

(^\ 

z+  \ 

Im 

= 

r  ,  z  = 

Im 

z-  J 

[o) 

[z-  ) 

(19) 


If  J,,,  denotes  the  m  xm  block  of  J  conformable  with  D^{k)  and  then 


T  =  D^{k)^J^,  +  D^{l). 


The  proof  is  so  similar  to  the  proof  of  Theorem  1  that  we  omit  it.  We  have 
allowed  for  the  fact  that  and  Z)_  need  not  have  the  same  number  of 
blocks. 

To  use  Theorem  6  to  approximate  an  eigenvector  suppose  that  J  is  nearly 
singular.  Compute  all  well  defined  T  and  find  one  with  a  minimal  singular 

o 

value.  Call  it  F-  Let 


r  V  =  U<T,ntn, 


ll«ll  =  llt^ll  =  1, 
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define  the  minimal  singular  triple 


(^mm  j  ^5  ^  )  of  F-  Then,  from  (19) 


J{Zv)  = 


0  \ 


oj 


^min* 


If  (^min  is  small  enough  then  Zv  is  &  good  initial  approximation  to  an  eigen¬ 
vector  of  J. 

It  is  not  hard  to  verify  that,  for  an  unreduced  even  order  J,  if  diag{J)  =  0 
then  diag{J~^)  =  0.  In  this  situation  a  block  factorization  with  2x2  blocks 
is  needed  to  ensure  that  J  =  L+D+U+  =  U-D-L-.  It  then  turns  out  that 


D+  =  D-  =  block  diag{J)  =  block  diag[J  ^ 

where  bdiag{M)  is  the  block  diagonal  part  of  M.  Now  the  set  of  F  matrices 
in  Theorem  6  may  give  no  guidance  for  computing  an  eigenvector.  That  is 
not  quite  true  because  we  may  infer  that  our  eigenvalue  approximation,  0,  is 
not  closer  to  one  eigenvalue  than  to  any  other  and  that  is  useful  information. 
In  fact  the  unreduced  J's  with  diag[J)  =  0  have  eigenvalues  in  i  pairs  and 
any  tiny  ±  pairs  may  be  found  efficiently  by  the  method  described  in  [13]. 
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5  5.5  6  6.5 

Figure  2:  Ritz  values  for  near  6 
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Figure  4:  log  7;  negligible  eigenvector  entries  from  42  to  105 


